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Abstract. Gaussian Quadrature is a well known technique for numerical in- 
tegration. Recently Gaussian quadrature with respect to discrete measures 
corresponding to finite sums have found some new interest. In this paper we 
apply these ideas to infinite sums in general and give an explicit construction 
for the weights and abscissae of Gaussian summation formulas. The abscissae 
of the Gaussian summation have a very interesting asymptotic distribution 
function with a (cusp) singularity. We apply the Gaussian summation tech- 
nique to two problems which have been discussed in the literature. We find 
that the Gaussian summation has an extremely rapid convergence rate for the 
Hardy-Littlewood sum for a large range of parameters. For functions which 
are smooth but have a large scale a the error of Gaussian Summation shows 
exponential convergence as a function of summation points, n: the error de- 
creases like ~ nexp(— 4?i^/7ra). The Gaussian summation achieves a given 
accuracy with n ~ ^/a points whereas other summation schemes require at 
least 71 ~ a function evaluations. 



1. INTRODUCTION 

Scalar sums of the form 

oo 

(1) E 

n— — oo 

with n integer and / being a continuous function appear in nearly every context 
of mathematics and physics. In many problems, either because / is for example 
the solution to another complicated nonlinear problem or there exists no known 
analytical expression for the sum, one has to resort to numerical summation tech- 
niques. Because of its importance many techniques have been developed for a fast 
evaluation of Eq. [T] which depend on the asymptotic behavior the function f{x). 
A good summary of the classical techniques can be found in [25. Newer develop- 
ments include application of Levin-type convergence transformation schemes ([9J 
for a review see [H|). In the application of these schemes the sequence of partial 
sums Sn = X]m=i /('^) Calculated and then extrapolated to n ^ 00. However if 
the function / has a large intrinsic scale so that the asymptotic behavior is reached 
only at very large values of n then these schemes are inefficient because a large 
number of terms has to be accumulated before the asymptotic behavior is reached 
and not applicable because rounding errors prevent a reliable extrapolation. So the 
problem of evaluating the sum Eq. [T] efficiently and reliably for problems with a 
large intrinsic scale is interesting. 
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In this paper we discuss a summation scheme based on ideas related to Gauss 
Quadrature. The basic idea is to replace the sum Eq. [T]by another sum 

oo N 

(2) E /H«E^fc/(^fc) 

n— — oo k—1 

where the weights Wk and abscissae Xk are chosen in such a way as to approximate 
closely the sum Eq.[T]for a large class of functions for a relatively small number N. 
Because the values of Xk are not fixed to be integer numbers the hope is that the 
asymptotic behavior at large x is captured by the sum Eq. [21 The central question is 
what are the optimal weights, Wk and points x^? This is the question we are going 
to address in this paper. The plan of the paper is as follows: First we derive the 
moment generating Weyl function for sums of the form J2ven l/^^/(V^^^) use 
the Fade approximants of the Weyl functions to obtain the orthogonal polynomials 
related to these sums. For numerical application we give the Jacobi matrix for 
generating the weights and moments for Gaussian summation formulas and give 
two examples of the application of the Gaussian summation. Finally we derive the 
asymptotic distribution of the zeros of the above orthogonal polynomials. 

Similar schemes have been reviewed recently by Engblom [3]. Milovanovic and 
Cvectkovic [Tl] studied the convergence properties of Gaussian summation schemes 
for sums of the form ^^^qP~'^ f {aq^'^), with > 1. The earliest reference 

to Gaussian summation can be found in the book by Uvarov and Nikiforov as an 
example using Chebyshev polynomials for finite sums. 

2. MOMENT GENERATING FUNCTIONS AND RECURSION RELATION 

We define a linear functional L[p], acting on the space of polynomials V, using 
some discrete measure A. 



(3) L[p]^ J pdX 

The focus of the paper is on sums of the type 

where fi is a subset of Z. The moments are defined as 

uen 

The moment generating Weyl function for Eq. 4 is given by 



z — X ^ — ' z V z ^ — ' 

i^en fc=o 

Weyl function can be calculated expHcitly for ft = Z\{0}: 

(7) a>(.)^i-^cot(^^ 

Eq. [7] serves as starting point for the deriving the orthogonal polynomials for the 
corresponding weight function. The Fade approximant of the moment generating 
function are closely related to the orthogonal polynomials of the corresponding 
weight function [l] . Here we derive an explicit expression for the Fade approximant 
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of Eqs. [7| which is equivalent to determining the orthogonal polynomials for the 
weight function Using the continued fraction expansion of cot(j;) we can 
express the generating function as: 

2 



(8) ^{x = tt/Vz) = 1 - a::cot(a:) = 



1 



3 

1 ' 



^ a;V63 



1 



The numerator of the nth Pade approximant corresponding to these continued frac- 
tions, Pn{x), is denoted by Rn{x) and the denominator by Sn{x). The denominator 
and numerator obey the recursion relations j2]: 

(9) Rn+l{x) = Rn{x) + Cn+i X^ Rn-l{x) 

(10) Sn+i{x) = 5„(x) + c„+i 5„_i(a;) 

for n = 0,1,2... with the initial condition for the recursion are = 0, 

Ro{x) = cqx^, S-i{x) = 1 and 50(2;) = 1. The coefficients c„ are given by cq — 1/3, 
c„ = -l/(2n + l){2n + 3), 71 = 1, 2, . . .. 

Surprisingly enough the recursion relation for Rn{x) and Sn(x) can be solved 
analytically. Let us discuss the recursion for Eq. [H First note that the recursion 
for the Bessel functions, of half-integer index n -f 1/2 are given by: 

(11) Zn+i/2{x) - '^^^'^ Zn+l/2{x) + Zn-l/2{x) = 0. 

Splitting of the asymptotic behavior for large n we define a„ (x) by requiring that 

/2\" 1 

(12) Z^+,/2{x)^i-\T{n+-)ar.{x). 
The coefficients a„(a;) obey the recursion relation: 



(13) 



An^ - 1 



This is precisely the recursion relation for the Sn{x) if we shift the index n by 2. 
Thus the recursion relation for Sn{x) has a solution of the form: 

rjf.n+2 

(14) S,,{x) = 5^ [A{x)J^+^/2{x) + B{x)Y^+^/2{x),] 



where A{x) and B{x) are two unknown functions which only depend on x and 
can be determined from the initial conditions. The orthogonal polynomials with 
respect to the functional Eq. |4]are given by denominators of the [M — 1/Af] Pade 
approximant [1] which are {iS'_i(a;), <S'i(a;), S^{x), S5{x), . . .}. Transforming back to 
the original variable z = {n/x)"^ we define the polynomials 

(15) s„(z) = Z"52„-1 

These are the sought after orthogonal polynomials for the functional, Eq. [H Spe- 
cializing to the initial conditions an explicit solution for Snix) can be obtained: 

(16) Sn{x = vr/Vz) = ~ p^^^ 5^2"+5/2 [cos(3;)'^n+5/2(a;) +sin(a;)y„+5/2(a;).] 
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This result will be useful to determine the distribution of zeros explicitly. The first 
few polynomials for n = 0, 1, . . . are given by 



So{x) 


= 1 






Si{x) 


= 1 - 


1 2 
X 

15 




S2ix) 


= 1 - 


2 2 

— x 
21 






= 1 - 


1 2 


1 4 
X 

945 


Siix) 


= 1 - 


±x^ 
33 


495 



We note that using Eq. [16] it can be shown easily that the polynomials s„ (x) are 
related to Bessel polynomials, yn{x), of imaginary argument [Hj: 



(17) 



s„(x) 



:Re ( y2«+i ( 



2 4"r(2n-, 2, 

Below we will derive simple three-point recursion relations for the polynomials 

Sn{x). 

3. Calculation of the weights and abscissae 

The orthogonal polynomials with respect to the functional Eq.[3]are {5*0, Si, S3, S5, 
It is useful to obtain a recursion relation for the odd numbered polynomials only. 
This can be achieved by eliminating the even numbered polynomials from the re- 
cursion. After some algebra we find: 
(18) 



Sn+3{x) = 1 



2x^ 



Sn+lix) 



-Sn-lix). 



(2n + 5)(2n + 9)y (2n + 3)(2n + 5)2(2n + 7) " 

This can be brought into the usual three-term recursion relation, Sk+i{x) = {x — 
afc)sfc(x) — 6fcSfc_i(x). The coefficients of the recursion for fc = 1, 2, . . . are given by: 

. , 2^2 
(19) 

(20) 



ak 



bk. = 



(4fc + 1)(4A:-H5) 

^4 



i4k-l){4k + lf{4k + 3)' 

The first coefficients have to be determined from the initial conditions and are 
Go = and bo is arbitrary. With these coefficients the problem of finding the 

weights and abscissae of the Gaussian summation formulas reduces to the standard 
procedure used to determine the analogous quantities in Gaussian integration [6]. 
The Jacobi matrix for determining the weights and abscissae is given by 

ao Vh 

Vh 0,2 Vbs 



(21) 



Jn — 



\/bN-i aN-1 
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It is easy to calculate the weights and abscissae using an implementation of the 
Golub and Welsh algorithm [6], either "gauss. m" from the QCP collection of Gautschi 
[4j, or the ^awco/ routine from the Numerical Recipes [T^. In terms of the weights 
Wn and abscissae a;„ determined from the Jacobi matrix, Eq. [2T]the Gaussian sum- 
mation procedure is given by: 



(22) = ^^,/(^,) + ,„ 

j/GO ^ ^ fc=l 

with the error estimate e„ (see e.g. [3J). The 

(23) e„ = K,.ig^ 

where ^ G C{n) and (7(0) is defined as the smallest connected set containing ft. 
The constant Kn can be calculated using Eq. (2.17) in [3]: 

(24) Kn = ||P„||' = ^(4n + 3)7r4"+3l6-("+i) T'^ (^^ + ^ 

Although the error estimate Eq. [23] gives an upper bound it is not very useful to 
estimate the convergence of the Gaussian summation formulas. For the specific 
examples below we can sharpen the error estimates and achieve a more general 
understanding of the convergence of the Gaussian summation formulas. 

4. Two EXAMPLES 

We apply the result of the preceeding section to two examples. The first is the 
Hardy-Littlewood function which has interesting properties in itself and is defined 



(25) H{x) = £ 



sin(a::/fc) 



k 

k=l 

This summand has an expansion in terms of 1/fc^ for /c/a; 1. For large values 
of X the summation is not expected to converge before the k x. The numerical 
evaluation of this function has been discussed by Gautschi [H \E\ . The methods 
considered in these references are Gaussian integration of the Laplace transforma- 
tion of the sum and the application of Euler-MacLaurin asymptotic summation the 
ehminate the first asymptotic terms 0(l/fc^) and 0(1/ k'^). 

We have applied our method to problem Eq. [25]for x = 100. The summation 
converges rapidly as shown below and for n > 14 the error in the sum is dominated 
by the error in the weights and points. One observes that the differences of the 
value of the Gaussian summation for different n is fiuctuating around the machine 
e. We find that the Gaussian summation scheme achieves a high accuracy using 
very few function evaluation of the summand. Basically the accuracy is Hmited 
by the accuracy of the Gaussian weights and points used in the evaluation of the 
summation. For example using fifteen terms the Gaussian summation gives the 
result for the Hardy-Littlewood sum with an error less then lO"^"* for all x G 
(0, 100). To compare specifically to the results by Gautschi in reference [4], Table 
3.21, with the largest value of x, x = 40, we compare the number of evaluations 
needed to achieve a relative accuracy of eo = 2.22 x 10~ as stated there. According 
to this reference the smallest number of points in the integration of the Laplace 
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transformed sum for this accuracy is iV = 39 so a total of 39 x 70 (for the integration) 
function evaluations are needed. With the Gaussian summation method a total 
n = 8 function evaluations are needed to achieve the same relative accuracy for 
H(AO). We summarize our results in Table [TJ The symbol — in the table indicates 
when the relative error in the summation was less then 10~^^. The error is plotted 
in Fig. [T] for X = 100 as a function of n. Gaussian summation give results already 
at a relatively small number of terms which does not grow like the scale where the 
function reaches the asymptotic behavior. 

Next we derive an analytic expression for the error estimate confirming the em- 
pirical findings. Obviously 



(26) j:4f] = Y,^kfixk)^ 

fc=i 



1 

27ri 



1 - 



-R2n-l(7I'C) 



1 



/ 



1 



where C is a contour enclosing the zeros of iS'„(7rC). The error estimate A„[/] = 
Sn+i[/] ~ can be simplified using the recursion Eq. [9] 



(27) 
with 

(28) 



A„[/] 



1 

2Td 



^ 27r(4n + 3)^4«+i(.4«+i 



1 



24n+5r2(2„ + 5/2) ^an+il^O-San-iW) 

The asymptotic behavior of the error is determined by Eq. [271 for n ^ 1. Using the 
asymptotic form of the Bessel function on the imaginary axis [lO^ and introducing 
J/ = 2n + 5/2 we obtain as the leading contribution: 

(29) A„(zi/r) « S^/e*'^^''^) 
with 

(30) $.(z) 
and 

For the Hardy-Littlewood function the error can be estimated as: 



E 




(^.±1)77 



■q{T) = v/t2 + 1 + log i^- 



(32) 



4i/ f 

— / dC exp($^(C)) sinh 
71" Jc 



This integral is completely dominated by the saddle points on the imaginary axis 
which appear for 2iP /ira ^ 1 and are located located at 



(33) 



±- 



Ve^iK 2(e-2)z. 



o 



1 



with the direction of steepest descent being perpendicular to the imaginary axis 
and ^ is given by ^ = 2v^ /ira. The evaluation of the integral at the saddle point 
yields: 



(34) 



2tt 



$"(i^Tmin) 



cxp ($^ (i^Tmin)) sinh 
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For 2v'^/TTa ^ 1 this expression can be simplified further to finally yield the as- 
ymptotic error estimate: 

(35) A„«2V^e--'^/'' (l^ 

We have plotted this error estimate in Fig. [TJ The agreement of the saddle point 
approximation Eq. [34] and the analytic asymptotic expression Eq. [35] with the nu- 
merical result is quite reasonable. This expression confirms that for smooth func- 
tions the Gaussian summation has a very rapid rate of convergence. The rate of 
convergence is determined by the scale of the function - in this case a at which 
the asymptotic behavior sets in. However in the Gaussian summation the rate of 
convergence is determined by the ratio a/v"^ so that convergence already sets in 
aX V K ^/a. This is a more rapid convergence rate than any of the extrapolation 
schemes for sums proposed in [5j. 

The second example serves to illustrate this point more clearly. We consider a 
sum known explicitly [10]: 

oo ^ 

(36) G(a)= ^2^ = ^^otM™), 

fc— — oo 

with a £ M>o- We have appHed Gaussian summation to Eq. [36] Again we observe 
rapid convergence even for values of a as large as 1000. The numerical convergence 
pattern is completely regular. In this case the expression for the error estimator 
A„ can be given without approximation because the contour integral can be done 
exactly by deforming the contour and is given by the residues of the integrand at 
zLia: 

(37) A„.-L/,CA„(C)^.= M^ 

Zm J (J a'^ -\- Q'^ a 

For large values of a using the asymptotic expression for An for v'^ ^ tto and v 
this expression simplifies to: 

/ 1/2- 

(38) A„ « 8i/ exp 

\ 7ra 

We contrast this behavior with other well known techniques for acceleration of 
sums like Richardson extrapolation [5]. Using the Euler-MacLaurin sum formula it 
can be shown that the partial sums behave asymptotically as: 
(39) 

^ , , 1 ^ 1 7rcoth(7ra) 2 1 - 1 60" - lOa^ + 1 
G„(a) = ^+2E ^TTF = -^-n^V-^^^-V^ +' 

k—1 

The A^'th Richardson extrapolation for the partial sums at some n is given by: 
(40) TIn[G\ - 2^ kHN-kV ^n+k{a). 

k=0 '' 

The asymptotic behavior of the sum, Eq.[36l is only reached when \k\ « \a\ and 
only then convergence of the partials sums starts. Applying Richardson extrapola- 
tion the terms in 1/n are successively removed. However because the coefficients 
of the n"™ terms are approaching a™ the convergence is very slow and the as- 
ymptotic behavior only sets in when n a. This can be seen in the plot of the 
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error in the Richardson extrapolations of Eq. [36l Fig. [31 The relative error is given 
by ATZN[G]{a) — \{TZN[G](a) — G{a))/G{a)\. It is clear that Gaussian summation 
achieve the same accuracy with a substantially smaller number of function evalua- 
tions. In fact the Richardson extrapolation shows spurious convergence for partials 
sums with N ^ a. 



5. Distribution of zeros of the orthogonal polynomials 

The distribution of zeros of the denominators of the Fade approximants Sn{x) 
turns out to be quite unusual. According to Eq. [I6]the zeros of Sn{x) for a given 
odd integer number N are determined by the equation: 

(41) cos(a::) J„+5/2(x) + sin(a;)y„+5/2(x) = 0. 

We are interested in the asymptotic behavior of the zeros at large values of n. Let 
a: > 0. The function is symmetric with respect to a; —> —a; so for each zero xq there 
is a corresponding zero at —Xq. We define v = n + 5/2 and t = x/v. The condition 
for a zero of the denominator apart from non vanishing prefactors is: 

(42) cos{vt)J^{vt) + sin{vT)Y^{uT) = 0; 

In the limit j/ ^ oo we have to differentiate between two regimes, r < 1 and r > 1. 
For the first case, r < 1 the leading asymptotic behavior of the Bessel function 
corresponding to a; < is given by [101 [7]: 

1 exp(i/(%/l - r2 - arcosh(l/r))) 

(4d) Ju{VT) = 



(l-r2)i/4 



and 

(44) Y,{vt) 



-2 exp(— j/(-\/l — — arcosh(l/T))) 

^ (1-t2)1/4 



The argument of the exponent of the second term is negative in the interval (0, 1). 
The first term is exponentially suppressed compared to the first term. The zeros of 
the denominator for t G (0, 1) are therefore determined by the expression: 

(45) sin{vT) = 0. 

From this equation we conclude that vt = ttto with m integer of the same order as 
V so it is useful to introduce a = m/u, so that in this regime r = ttct. 

Next we consider the case r > 1. The leading term in the asymptotic expansion 
of the Bessel functions takes a different form: 

/ 2 

(46) Jyiyr) = W — cos(z/(v/r2 - 1 - arccos(l/T))) 

V 'KV 

/ 2 

(47) Y^{vt) = \ — sin(i/( - 1 - arccos(l/r))). 

V TTV 

For this case the leading asymptotic behavior of the zeros of the denominator have 
to be determined by the following equation: 



(48) cos(z^(t — VT^^-l + arccos(l/r)) = 0. 
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In this regime the zeros are solutions of t — y/r^ — 1 + arccos(l/r) = na. To 
summarize the asymptotic behavior of the zeros is given by: 



Note that a approaches 1/2 in the Hmit r ^ oo. This is the expected because we 
have considered the case a; > only. In fact the distribution of zeros is symmetric 
in the range a e (—1/2, 1/2) since the Sn{x) depend x'^ so only half of the zeros lies 
above 0. We enumerate the positive zeros by Xk, k £ {1,2, . . .n}. The n negative 
zeros are given by —Xk- The asymptotic behavior of the zeros is surprising: For 
large n the first [n/n] zeros are with high accuracy given by x = tt, 2tt, Stt, .... How- 
ever for k > [n/Tr] the distance between two consecutive zeros become larger and 
eventually for very large Xk becomes very large and behaves like Xk ''-^ j-Kiy — k). 
The asymptotic behavior has been verified in numerical determination of the zeros. 
This result is unusual because a singularity appears in the distribution of zeros 
in the range of integration which is not present in the moment generating func- 
tion. The Gaussian summation achieves its accuracy by evaluating [n/ir] points at 
the summation points and using the rest of the points to explore the asymptotic 
behavior of the function which is being summed. 



The general problem of summation of smooth functions is considered. Summa- 
tion formulas in analogy to Gaussian quadrature have been derived. The Pade 
approximants of the moment generating functions have been used to derive the 
expHcit solution of the recursion relation. We have used this to give an expHcit for- 
mula for the corresponding orthogonal polynomials. The asymptotic distribution 
of the zeros of these polynomials is derived. We have shown that the distribution of 
zeros shows an unusual behavior exhibiting a cusp. We demonstrate that Gaussian 
summation needs a substantially smaller number {^/n compared to n) of function 
evaluations for smooth functions with a large scale than other common techniques. 
In particular higher order Richardson extrapolation will lead to rounding errors 
spoiling the convergence. 

It would be interesting to investigate Gaussian summation for different kind of 
linear functionals like X)iyGZ\{o} 1/^"^ fi^/^'^)- Summation schemes for summation 
over odd numbers (Matsubara frequencies) relevant for semiclassical field theories 
in physics will be discussed in a separate paper |12j . Another interesting point is the 
development of Kronrod schemes for summation which we are currently working 
on. 
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(49) 




r < 1 



T > 1. 



6. Conclusions 



7. Acknowledgments 
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n 


AH{1) 


AH{5) 


\ // M 


AH{20) 


AH{40) 


Aff(lOO) 


2 


8.73 • 10"^ 


1.58 • IQ-^ 


1.9 • 10° 


4.47 ■ 10-1 


5.06 • 10-1 


5.61 • 10-1 


3 


- 


2.53- IQ-*^ 


1.02 • 10-2 


9.55- 10-1 


2.29 • 10-1 


3.29- 10° 


4 


- 


3.66- IQ-" 


3.3 • 10-6 


1.67- 10-2 


1.28 • 10° 


3.29- 10° 


5 






1.47- 10-1° 


2.29 • lO-'"" 


2.48 • 10-1 


4.51 • 10-1 


6 








5.19 • 10 


o ri A in — 3 

3.04 • 10 


2.53 • 10 


7 








2.85-10-13 


5.89-10-6 


2.77- 10° 


8 










2.8- 10-^ 


1.09- 10° 


9 










4.19- 10-13 


4.46- 10-2 


10 












3.87- 10-^ 


11 












1.02- 10-s 


12 












1.01- 10-9 


13 












4.07- 10-13 


14 












2.51 • 10-" 


15 















Table 1. Results for the relative error in evaluation of the Hardy- 
Littlewood sum. 
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10° 

10-= 

10-^° 
^" 10- 
10-2° 
10-2= 
10-3° 

2 4 6 8 10 12 14 16 18 
n 

Figure 1. Relative error in the evaluation of the Hardy- 
Littlewood sum for x = 100. The squares (□) are indicating the 
numerical error in the Gaussian summation, the crosses (x) are 
indicating the saddle point approximation of the error integral and 
the solid line gives the asymptotic behavior of the former for large 
n. The dashed line gives the machine e for comparison. 
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Figure 2. Relative error of the Gaussian summations of the ex- 
ample 2, sum Eq. [SH for a = 1000. The solid line is the error 
estimate given in the text 8n exp(— 4n^/7ra)) and the dotted 
line is giving the machine e. 
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Figure 3. Relative error of the Richardson extrapolations of the 
example 2, sum Eq. [36l for a = 1000. Note the appearance of 
rounding errors in A7?,4 [G] (n) . 
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Figure 4. Asymptotic density of zeros ( ii{a) = ^(cr)). The dots 
give the approximate density of zeros l/(a;„+i — a;„) for N = 32 
plotted as a function of cr = n/N. The dotted line indicates a = 
I/tt. 



